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We investigate the dynamics of Rydberg electrons excited from the ground state of ultracold 
atoms trapped in an optical lattice. We first consider a lattice comprising an array of double-well 
potentials, where each double well is occupied by two ultracold atoms. We demonstrate the existence 
of molecular states with equilibrium distances of the order of experimentally attainable inter-well 
spacings and binding energies of the order of 10 3 GHz. We also consider the situation whereby 
ground-state atoms trapped in an optical lattice are collectively excited to Rydberg levels, such that 
the charge-density distributions of neighboring atoms overlap. We compute the hopping rate and 
interaction matrix elements between highly-excited electrons separated by distances comparable to 
typical lattice spacings. Such systems have tunable interaction parameters and a temperature ~ 10 4 
times smaller than the Fermi temperature, making them potentially attractive for the study and 
simulation of strongly correlated electronic systems. 

PACS numbers: 37.10.Jk, 34.20.Cf, 32.80.Ee 



I. INTRODUCTION 

Recent advances in the trapping and manipulation of 
ultracold atomic gases have provided experimentalists 
with the ability to coherently control large numbers of 
atoms. Two areas that have become the focus of ex- 
perimental efforts of late are the use of ultracold atoms 
for the formation and manipulation of molecules [l], Q 
and the creation and manipulation of Rydberg atoms 
in optical lattices [3J. In this paper we study whether 
combining these areas might lead to the production of 
diatomic molecules whose nuclear position is fixed by an 
optical lattice (see Figs.QJ,-b). In particular, we examine 
the properties of ultralarge dimers with equilibrium dis- 
tances of the order of typical lattice spacings and binding 
energies of the order of 10 3 GHz. We also investigate the 
prospect of using systems of interacting Rydberg atoms 
to simulate Fermi systems. 

Molecules have a far richer energy structure than 
atoms, and can also have stronger long-range interac- 
tions, a feature which offers new possibilities for quantum 
control (see e.g. Q). However, cooling molecules is no- 
tably difficult, since the absence of closed electronic tran- 
sitions prevents the use of standard laser cooling proce- 
dures. An attractive approach to producing translation- 
ally cold molecules is thus to form them from pre-cooled 
atoms by way of photoassociation or magnetic resonance 
techniques. In recent years several classes of ultracold 
molecules have been predicted and produced, includ- 
ing Feshbach molecules, Efimov trimers, and the famous 
'trilobite' molecules, which are composed of a highly ex- 
cited Rydberg atom interacting with another atom in its 
ground state [3, @ . The existence of long-range molecules 
composed of two Rydberg atoms stabilized via dipole- 
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FIG. 1: (color online) An optical lattice with pairs of atoms 
well separated from each other is initialized (a) . A laser pulse 
transfers each pair of atoms to a molecular state with a very 
large internuclear distance, (b) Density plot of a typical di- 
atomic molecular wavefunction on the x-z plane (y = 0) 
in the relative coordinates of two electrons in highly ex- 
cited np z states, (c) The outer electron of each ground-state 
atom trapped in an optical lattice is transferred into a Ryd- 
berg state, such that the charge-density distributions between 
neighboring atoms overlap. The electron hopping rate t be- 
comes non-zero, and the interactions between electrons can 
be described by the parameters U (on-site interaction), V, W 
and X (off-site interactions). 



dipole interactions has also been predicted 0, H, H|. 
These ultracold long-range Rydberg molecules have bind- 
ing energies of the order of several hundred MHz (a factor 
of ~ 10 6 greater than the typical temperature of ultra- 
cold atoms, but weak by molecular standards), a lifetime 
expected to be similar to that of Rydberg atoms, and 
equilibrium distances of the order of 10 cto (with ao the 
Bohr radius) [8| . Recent advances in high-resolution mi- 
croscopy techniques may offer new possibilities to study 
the spatial structure of such molecules [10(. Theoretical 
treatments of Rydberg molecules to date have primarily 
dealt with the regime in which the internuclear separa- 
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tion is greater than the Le Roy radius. In this regime the 
overlap between the atomic charge-density distributions 
is vanishingly small 0, @|- In contrast to these works 
we consider the regime where the internuclear distance 
of the molecular dimer is comparable to typical lattice 
spacings, but smaller than the Le Roy radius, in which 
case the overlap of the charge distributions and effects 
such as the exchange interaction must be taken into ac- 
count. 

The simulation of condensed-matter systems is another 
area to which optical lattices are uniquely suited. Opti- 
cal lattices have extremely flexible geometries, and by 
a suitable choice of laser configuration, any desired lat- 
tice structure can be created [111, [HJ . Fermions loaded 
into an optical lattice can be used as a model of elec- 
trons in a solid, with the significant advantage of having 
a periodic potential that is defectless and fully customiz- 
able [131 ]. However, the ratio between the currently at- 
tainable temperature T of fermionic atoms trapped in 
an optical lattice and the Fermi temperature Tp of the 
system is T/Tp ~ 0.25, preventing the clean observation 
of such interesting phenomena as the BCS transition or 
the emergence of certain types of anti-ferromagnetic or- 
der [l3|, [lj, . In the final section of this paper we con- 
sider the possibility of collectively transferring a popula- 
tion of ground-state atoms trapped in an optical lattice 
to a highly excited Rydberg level, such that the charge- 
density distributions of neighboring atoms overlap (see 
Fig. [TJ:). We extend the model developed in the first 
section to compute the hopping rate and interaction pa- 
rameters between the highly excited electrons, and deter- 
mine the dependence of these quantities on the initial lat- 
tice spacing. While the implementation of these systems 
is experimentally challenging and beyond the scope of 
this paper (l6| . we note that experimental temperatures 
far smaller than the Fermi temperature {T/Tp ~ 1CP 4 ) 
might be readily attainable, making the realisation of 
such systems a promising new approach to the quantum 
simulation of interacting fermions. 

The paper is organized as follows. In Sec.|TI]we present 
the model that we subsequently use to investigate the 
existence of attractive molecular potentials with equilib- 
rium distances below the Le Roy radius. We approximate 
the lifetime of these molecules and calculate their equi- 
librium internuclear distance for different values of the 
principal quantum number of the electrons. In Sec. IIIII 
we use the model introduced in the first section to com- 
pute the hopping rate and interaction parameters of the 
Rydberg electrons and discuss the possibility of using 
these setups to simulate condensed-matter systems. We 
conclude in Sec. [TV] 



II. MOLECULAR POTENTIALS OF 
ULTRA-LARGE RYDBERG DIMERS 

In this section we present a method of evaluating the 
energy of highly excited diatomic molecular states whose 



equilibrium distance R e is below the Le Roy radius. We 
develop an efficient method of computing one- and two- 
center molecular integrals over electronic orbitals with 
large principal quantum numbers n and calculate the 
depth and equilibrium position of molecular potentials 
accessible from the ground state via photoassociation for 
values of n up to n = 35. We estimate the radiative 
lifetime of the highly excited molecular potentials. 

Quantum chemistry methods commonly used within 
the Le Roy radius typically become computationally ex- 
pensive when dealing with Rydberg molecules due to the 
size of the basis set required to describe the spatially dif- 
fuse Rydberg orbitals. In order to obtain qualitative re- 
sults in this computationally challenging regime we take 
a simple wavefunction ansatz inspired by the Heitler- 
London treatment of the hydrogen molecule. Due to the 
simplistic form of the molecular wavefunctions used we 
do not expect the method presented to produce quan- 
titatively accurate results; rather, our aim is to obtain 
qualitatively correct results and insights within a regime 
where common methods of calculating molecular poten- 
tials fail or become computationally costly. 



A. Model 

We consider an optical lattice potential forming an ar- 
ray of double-wells separated from each other by suffi- 
ciently high optical barriers that each double well can 
be considered as an isolated system [l7|. We assume 
that each double well initially contains two 87 Rb atoms 
in their ground state, and that the potential is suffi- 
ciently deep that the Wannier function associated with 
each atom is well localized in one half of the well. Because 
alkali atoms have only a single valence electron, their en- 
ergy levels are described by the same quantum numbers 
as those of the hydrogen atom. If the ground-state atoms 
are excited to Rydberg levels with n > 30 by applying 
a laser pulse to the system, the charge distributions of 
atoms in the same double well will overlap, and under 
certain conditions a stable molecular state will be formed. 
We aim to study the basic properties of these molecules 
and explore whether they may be produced in an optical 
lattice with experimentally realistic parameters. Since 
the temperature of the system is very low, the velocity 
of the electrons - even in highly excited states - is much 
greater than that of the trapped ions, and the electrons 
respond almost instantaneously to displacements of the 
ions [8[. Consequently, we will calculate molecular po- 
tentials in the Born-Oppcnhcimcr (BO) approximation 

nam. 

Neglecting interactions between atoms belonging to 
different double wells, the energy of a homonuclear 
molecule formed by the atoms in a double well is given 
within the BO approximation by the eigenvalues of the 
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Hamiltonian 

\R ri2 n A n B r 2 A r 2 B J 

where r; (i — 1,2) is the position of the i th electron; 

(£ = A, B) is the distance between the i th electron 
and atomic center £; ri 2 is the distance between the two 
electrons; jo = e /(47T£o); where e is the permittivity 
of free space and e the elementary unit of charge; and R 
is the distance between the nuclei. In the BO approx- 
imation R is treated as a classical variable upon which 
the eigenvalues and eigenfunctions of the Hamiltonian 
depend parametrically. Following a similar approach to 
that of Ref. @, [2(|, we approximate the energy of the 
molecule by diagonalizing the Hamiltonian (TTJ) using as 
a basis the asymptotic electronic states 

<q' =-^ q <[A q (ri)B q ,(r 2 ) ± A q (r 2 )B q ,(n)], (2) 

where A p (ri) (B p (ri)) is a hydrogen-like wavefunction 
with quantum numbers q = (n, £, m) centered on nucleus 
A(B), and jV± q , = [2(1 ± |5 q , q - | 2 )]~ 1/2 is a normaliza- 
tion factor, where S q , q ' = J dr A*(r)B q '(r) is the over- 
lap integral [l|j]. The basis elements that are symmetric 
or anti-symmetric in the coordinates have an associated 
anti-symmetric (singlet) or symmetric (triplet) function 
of spins respectively. Each basis state (J2|) describes the 
two electrons being exchanged between the two nuclei, 
and corresponds to an ansatz for a valence-bond wave- 
function [13|. The diagonal elements of the Hamilto- 
nian |T]) using the states ([2]) as a basis correspond to the 
molecular potentials obtained using the Heitler-London 
method [IE EH. 

The hydrogenic wavefunction of the valence electron 
centered on atom £ located at is well approximated by 
the (unnormalizcd) function 

£ q (r)=P,(cos % )!^ (y"'e-/»- l fb t ^ w ». 

a V / k=0 

(3) 

Here = r^/ao with = |r — r^| the radial distance 
from atomic center £, 9^ is the angle between the vector 
and the internuclear axis, <f> is the azimuthal angle, 
P(.(x) is a Legendre polynomial, and n* — n — 5i is the 
effective principal quantum number, with Si a quantum 
defect whose value depends on the angular momentum 
quantum number £ and the atomic species (22[. The co- 
efficients b fe = b fc _i (n*/2k) [£(£ + 1) - (n* - k)(n* - k + 1)] 
are defined recursively with bo — 1 and fc ma v- is an integer 
satisfying n* - 1 - 1 < fc max < n* - £ [H E|] . The orbital 
defined in Eq. ([3]) is known as the asymptotic form of the 
quantum defect wavefunction; for 5i = 0, it is identical 
to the hydrogenic wavefunction, while for Si ^ it pro- 
vides an accurate description of a Rydberg electron in the 
mid- and long-range. It differs significantly from the ex- 
act quantum defect wavefunction only at short distances 



from the atomic cores, which is of little consequence since 
the interactions we consider here depend mainly on the 
outer part of the atomic wavefunctions. 

For n > 20, the overlap between the charge-density dis- 
tributions of atoms separated by distances smaller than 
the Le Roy radius (for q = q') is typically of the order 
of S q , q ' ~ 10 _1 — 10~ 2 . The charge-density overlap is 
therefore not negligible in this regime, and the approxi- 
mation of the Coulomb potential using the multipole ex- 
pansion [25L l2rl | is not suitable. The multipole expansion 
diverges quite significantly from the Coulomb potential 
even in the presence of small charge-density overlap, and 
so the results obtained using this approximation below 
the Le Roy radius are very uncertain [1, E3 • 

We consequently estimate the energy of the dimcr 
by evaluating the expectation value of the Hamilto- 
nian ([1]) in the states ([2]), a task which requires the 
evaluation of a number of integrals over atomic orbitals. 
These integrals fall into two classes: the one-center 
integrals, comprising the overlap integral S'q.q/, the 
Coulomb integral J q , q - = j f drB*(r)(l/rA)B q /(r) = 
jo / dr A q (r)(l/re)A q / (r), and the charge overlap inte- 

S ral K Lq> = jo I drA q (r)(l/r|)B q ,(r); and the Oo- 
center integrals 

^ P q p q ' = [A q A q < I A P A P ,] W™, = [A q B q , I ApBp,] 

(4) 

= [A q A q , |BpBp-] X$$ = [A q A q , |A p Bp/] 

(5) 

where 

[a/3\jv]=j J J dridr 2 a*(ri)/3(ri)-5-7*(r2Mr 2 ). 

(6) 

The evaluation of one- and two-center molecular inte- 
grals poses a considerable challenge for large values of 
the principal quantum numbers. The results of such in- 
tegrals using direct numerical integration converge ex- 
tremely slowly, and the answers so produced can suffer 
from dramatic losses of accuracy (a p henomenon known 
as numerical erosion; see e.g. [28j, 29]). Following an ap- 
proach suggested by M. P. Barnett, we use a computer 
algebra-based method to generate analytical formulae for 
the molecular integrals [23, E(J HH . We have found this 
approach advantageous for three reasons: (i) it permits 
the fast evaluation of molecular integrals with high prin- 
cipal quantum number to arbitrary accuracy; (ii) once 
the analytical form has been found, the evaluation of an 
integral for different values of the internuclear distance is 
instantaneous; (iii) the analytical expressions of the in- 
tegrals can be stored and re-used at little computational 
cost. 

However, even with the use of symbolic calculations, 
the evaluation of molecular integrals for large values of 
n is computationally demanding. In order to make these 
calculations tractable for large n we restrict the value of 
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Rule Functional Replacement 

Rl f~dxe- ax x k - (l/a 1+k )T 1+k (a) 

R2 \f^dxe^x k -> [(-(3)- k /0\\r 1+k (-l3)-r 1+k <J3)] 

TABLE I: Replacement rules used to generate analytical for- 
mulas for one-center integrals; I\ (x) is the incomplete gamma 
function with k £ N+. 



the projection of the electronic angular momentum along 
the internuclear axis to m = 0. This limits the range of 
molecular states that can be investigated to those of sym- 
metry and 3 S+; these are associated with the >F q ~ , 
and *f?~ , basis states respectively. In this way we have 
been able to compute molecular integrals involving wave- 
functions with principal quantum numbers up to n = 35. 

In order to compute the molecular integrals it is conve- 
nient to express the electron coordinates in elliptical coor- 
dinates (A, fi, </>) through the relations = (i?/2)(A ± /i) 
and cos0£ = [(1 ± A/z)/(A ± /J,)], where plus and minus 
signs apply to £ = A and £ = B respectively. The vol- 
ume element is given by dr = (R/2)rArsd\ d/i dip, where 
< (j) < 2ir, — 1 < fjb < 1 and 1 < A < oo. After rounding 
up the powers of in Eq. ([3]) to the next integer value 
(which only significantly affects the shape of the wave- 
function near the core) we find that the integrands of all 
of the one-center integrals can be written in the form 

e -aX e /3ti ^ qtjX 1 ^ , where qij are coefficients associ- 
ated with a given integral. By applying the replacement 
rules defined in Table [J every one-center integral can be 
converted into an analytical expression with parametric 
dependence on R. 

The evaluation of symbolic expressions for two-center 
integrals proceeds similarly although many more replace- 
ment rules are required. The general approach for two- 
center integrals consists of expressing the Coulomb po- 
tential in the form of a sum of polynomials (such as 
the Legendre or Neumann expansion) before applying 
replacement rules on the terms resulting from the suc- 
cessive integrations over the coordinates of the first and 
second valence electron. Using this method, exact for- 
mulae are obtained for Eqs. ([3]), but only approximate 
expressions may be found for Eqs. (O (see Appendix [A"| . 

B. Results 

Since we envisage producing the Rydberg dimers by 
applying a laser pulse to a system of ground-state atoms, 
we are primarily interested in the molecular states that 
are most strongly coupled to ground-state atoms by the 
dipole transition operator; namely, those of 3 £+ symme- 
try with a high p-character (£ = 1). We will therefore 
restrict our analysis to these states. However, for values 
of n > 10 molecular states of 3 E+ and symmetry 
become quasi-degenerate d, [26[ , and so the results pre- 
sented here apply to molecular states of either of these 
symmetries. 
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FIG. 2: (color online) Equilibrium distance R e (a) and poten- 
tial depth D e of the np + np states using the Heitler-London 
method (b). The dots denote the values calculated numeri- 
cally, while the curves correspond to the fitted functions men- 
tioned in the text. 



We have estimated the equilibrium distances R e and 
potential depths D e = E mo \(oo) — E mo \(R e ) of np + np 
molecular potentials of 3 E+ symmetry using the Heitler- 
London approach for values of the principal quantum 
number up to n = 35 (see Fig. [5]). We find that the 
equilibrium distance of these potentials follows the rela- 
tion R e ~ (1.628n 2 — 100.25) do, which is about a factor 
of four smaller than the Le Roy radius for all values of 
n (see Fig. [5^,) . For larger values of the principal quan- 
tum number (n > 35) we find that the equilibrium dis- 
tances predicted by our scaling law are comparable to 
(although ~ 20% smaller than) those found by Boisseau 
et al. in Ref . , who use the multipole expansion of the 
Coulomb potential to examine the same potentials below 
the Le Roy radius. 

The depth of the molecular potential for n = 35 is 
D e ~ 2000 GHz, approximately three orders of mag- 
nitude larger than the potential depth of ultra-large 
molecules with symmetry 1 n £l — 3 IL. bonded via dipole- 
dipole interactions studied in Ref. [8[ (see Fig. [2Jd) . We 
find that the molecular binding energy decreases expo- 
nentially with the value of the principal quantum num- 
ber according to the relation D e = exp(ai + a^n^ 2 + 
cv3?7r 3 )GHz where a% — —6.53, 012 = —24.44, = 
-20.51, fo = 0.14 and /3 3 = -223.45. Our predicted 
potential depths differ by up to an order of magnitude 
from those calculated using the multipole expansion; as 
Boisseau et al. point out in Ref. Q, this is probably a 
consequence of the inaccuracy of the multipole expan- 
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FIG. 3: (color online) Left axis: The 3 E U molecular potential 
(solid red line) associated with the wavefunction \& mo i(-R) re- 
sulting from the diagonalization of the Hamiltonian (see text) 
for n = 16. The dashed line shows the molecular potential 
for the np+(n— l)p bare state. Right axis: The contribution 
Cq,q'(R) of the q = (16, 1,0), q' = (15, 1,0) basis state to the 
molecular wavefunction * mo i(-R); A denotes the width of the 
region over which this contribution is dominant. 



sion below the Le Roy radius. Although the potentials 
produced using our method do not have the asymptotic 
R~ 5 behaviour expected from perturbation theory [26| . 
this may be corrected by carrying out a numerical in- 
tegration of the terms which were neglected during the 
symbolical calculations of the molecular integrals. This 
correction is feasible for small values of n, but very time 
consuming for larger values. However, a comparison with 
the results of exact numerical calculations of the same po- 
tentials up to n = 8 shows that our method reproduces 
the correct values of R e and D e within ~ 3%. 

Following the approach taken in Ref. 0], we diago- 
nalized the Hamiltonian using as a basis the molecular 
states ([2]) with a significant coupling to the np + n'p 
asymptote. For n' = n— 1 and n = 16 we used the states 
(n-l)s+(n + l)s, (n-l)s + ns, (n-l)s + (n + 2-k)d, 
(n - 2)d + (n - k)d, (n - 3)d + (n - k + l)d and 
(n - 1 - k')p + (n + k)p, with k,k' = 1,2. The di- 
agonalization procedure yields eigenstates of the form 
'I'moi(-R) = ^Cq, q '(iJ)>tq,q', allowing the determination 
of the contribution of a given molecular state ^q.q' to 
each eigenstate. Figure [3] shows the contribution of the 
asymptotic np + (n— l)p state to an eigenstate ^moi, and 
the potential curve -Emoi {R) corresponding to this eigen- 
state, for the case of n — 16. In a region of width A about 
R e the molecular potential is indistinguishable from that 
of the bare np + n'p state, and the associated wavefunc- 
tion has a very dominant np+n'p character (C qiq ' = 0.98 
with q = (16,1,0) and q' = (15,1,0); see Fig.©. This is 
due to the fact that the coupling between different states 
is generally weak (off-diagonal elements of the Hamilto- 
nian are typically ~ 10 _1 — 10~ 3 smaller than diagonal 
ones) and the strongest coupling is seen between states 
with identical angular momentum. 

We find that the size of the region A ~ 0.24n 2 ao 



is approximately a factor of three smaller than the 
width of the molecular potential associated with the bare 
np + (n — l)p states (approximately given by 0.73n 2 ao). 
This provides an indication of the nature of the wave- 
function of the molecular states accessible by the pho- 
toassociation of ground-state atoms initially trapped in 
an optical lattice. The distance adw between the two 
sides of a double well may be controlled by altering the 
laser parameters. The uncertainty in the initial posi- 
tion of the ground-state atoms is given by the width of 
their Wannier functions (see Fig. [XJ) , which is propor- 
tional to a = ad w /ln(V /En) 1 / 4 ], where V is the depth of 
the double- well and Er the recoil energy [33] ■ By setting 
&dw = Re and considering the scaling of the width of the 
np + (n — l)p molecular potentials, we find that 2a < A 
for lattice potential depths of the order of V = 30-40 Er. 
This suggests that using an optical lattice to enforce the 
initial position of the atoms before applying the photoas- 
sociation pulse might offer advantages in providing access 
to molecular states with a very dominant np + n'p char- 
acter. 

We have obtained an order-of-magnitude estimate of 
the lifetime of the photoassociated molecules by assum- 
ing for simplicity that the dominant decay mode is radia- 
tive dissociation into a pair of 87 Rb atoms. We neglect 
the rotational fine structure of the energy levels and any 
bound-bound decay channels and consider a transition 
between a bound state of vibrational quantum number v' 
and energy E v i and a free continuum state of wavenum- 
ber k" and energy Ek" — h 2 k" 2 / (2fi), where \i = m^b/2 
is the reduced mass of the molecule. The Einstein A 
coefficient is given by 
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(7) 

while the radiative lifetime of a single bound vibrational 
level is given by t v i = A~, , where 



A,jiyi dEi-ii . 



(8) 



Here h = 2irh is the Planck constant, c is the speed of 
light, ipl) h and ij/Q^ are vibrational wave functions for 
the discrete and continuum states respectively, and the 
dipolc moment D(R) is given by 



D{R) =-e dndr 2 + z 2 ] 



q/.q/ 



(9) 



where z\ and zi are the z-coordinates of the electrons 
and (\i and q/ denote the quantum numbers of the ini- 
tial and final states. For large internuclear distances the 
continuum wave function ipffi has the asymptotic form 



sm(k"R + i]) 



(10) 
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where r\ is the scattering phase shift, which at low en- 
ergies is given by r\ = —k"a s with a s the s-wave scat- 
tering length. In the initial molecular state we approxi- 
mate the true bound vibrational wavefunctions \jf§°(R) 
by the eigenstates of a harmonic approximation to the 
molecular potential centered about R = R e - Due to 
the asymmetric nature of the molecular potential curves, 
this approximation becomes progressively worse as v' in- 
creases. However, we expect that the ability to impose 
an interatomic spacing close to the equilibrium distance 
of the targeted molecular state before applying the pho- 
toassociation pulse will provide a measure of control over 
the range of vibrational levels occupied by the molecules, 
thus making the lowest vibrational levels the most rele- 
vant. 

Although there are a huge number of final states to 
which the molecule could decay, the dependence of the 
decay rate given by Eq. ((7]) on [E v i — Ey ) 3 favors transi- 
tions that involve the emission of a high-energy photon. 
We therefore consider only transitions that finish within 
the continuum above the 5s-5s potential curve, neglect- 
ing all other decay channels. The radiative lifetimes thus 
calculated for the ground vibrational state v 1 = are 
shown in Fig. 2J the lifetimes of higher vibrational lev- 
els (up to v 1 — 10) are the same to within ~ 5%. Also 
plotted for comparison purposes is the radiative lifetime 
r = rpn 3 of a free Rydberg atom, where To = 1.4 ns for 
87 Rb [33| . Our calculations indicate that for n > 17 the 
Rydberg dimers are more stable than the free atoms. For 
the highest molecular state considered (n = 35) we find 
lifetimes of the order of a few milliseconds, indicating 
that even if contributions from neglected decay channels 
were to reduce this lifetime by several orders of magni- 
tude, the system could still be successfully interrogated 
and characterised by short (nanosecond to picosecond) 
laser pulses. To avoid possible stimulated emission con- 
tributions to the radiative decay process, the laser fields 
forming the optical lattice could be switched off as the 
Rydberg excitation pulse is applied. The subsequent free 
expansion of the atoms would not limit the window of 
time within which the system may be characterized, since 
the atoms are expected to remain localized within a few 
lattice spacings for a time of the order of 10 /zs [34| . 
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FIG. 4: (color online) Radiative lifetimes for the v' = vibra- 
tional level in molecular excited states with principal quan- 
tum numbers n = 15-35 (black line). The radiative lifetime 
r = Ton 3 of a free Rydberg atom is plotted for comparison 
purposes (red line). 



that the band structure of such a system is well described 
by a tight-binding model, we will find that using the hop- 
ping rate calculated in the next section that the Fermi 
temperature of a half-filled lattice is Tp ~ 10 _3 K, some 
5 orders of magnitude larger than typical temperatures 
achieved in experiments with ultracold atoms. While the 
efficient population transfer of atoms to Rydberg states 
is experimentally very challenging and goes beyond the 
scope of this paper (see e.g. Ref. [lj| ), this figure im- 
plies that after the excitation step, the system tempera- 
ture could still remain well within the Fermi degeneracy 
regime (T <C Tp). This would enable the observation of 
interesting many-body phenomena currently inaccessible 
to other experimental setups involving fermionic atoms 
where the ratio T/T F » 0.25 (see e.g. Ref. 0, S3). In 
this section, we use the tools developed previously to 
evaluate the typical interaction parameters and the hop- 
ping rate of such a system, and investigate their depen- 
dence on the interatomic spacing fixed by the lattice. 



A. Model 



III. INTERACTIONS BETWEEN HIGHLY 
EXCITED ELECTRONS IN A LATTICE 

In this section we consider the situation where ground 
state atoms trapped in a regular optical lattice are col- 
lectively excited to a given Rydberg state such that the 
charge-density distributions of neighboring atoms over- 
lap. With overlapping charge-density distributions be- 
tween neighboring atoms, the valence electrons will tun- 
nel between lattice sites and interact with each other, 
mimicking the behaviour of electrons in metals. Such a 
system may have interesting applications for the quan- 
tum simulation of electrons in lattice systems. Assuming 



In this section we assume that ground-state atoms 
trapped in an optical lattice can be collectively excited 
to a given Rydberg state. We assume that the dynam- 
ics of the valence electrons in the lattice is restricted to 
only one spatial mode per site i that corresponds to an 
orbital <pi(r) = (r|0j) associated with a given Rydberg 
state localized at site i with two spin orientation. The 
Hamiltonian of the system is then given by 

i,j.cr i,j,k,l 
cr,cr 

(ii) 
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FIG. 5: (color online) Relative values of the interaction pa- 
rameters (a = V,t, X , W) of a Hubbard model of Rydberg 
electrons in a state with quantum numbers q = (30, 1, 0). (a) 
I = (b) I = 0, setting S q , q = 0; (c) £ = 1 (d) I = 1, setting 
5q, q = 0; (e) £ = 2 (f) £ = 2, setting 5 q , q = 0. From top to 
bottom we have V/U (green), t/U (blue), X jU (orange), and 
W/W (red). 



where c ia is a creation operator associated with the or- 
bital 4>i( r ) an d spin a, Uj describes the hopping of a par- 
ticles between sites i and j; and the inter-site interactions 
are given by 



V{i,j,k,l) = J J drdr>*(r)^(r')Kc(r-r')0 fc (r)0 ; (^), 

(12) 

where V ec (r — r') is the interelectronic potential. We as- 
sume here that next-nearest-neighbour hopping can be 
neglected, and consider only nearest-neighbour interac- 
tions. In this situation the hopping term t = t^j+i 
consists of the kinetic energy and Coulomb potential of 
neighboring ion cores 



terms of electronic orbitals is given by 

4>i = ipi ^ ■ 



(14) 



where ipi is a normalized electronic wavefunction with 
quantum numbers centered at site i, and Si is the overlap 
between the orbitals of site i and i + 1 [38| . In the regime 
where the terms proportional to Sf can be neglected, the 
functions (pi form an orthonormal basis set. By inserting 
the definition (Thi]) into Eqs. (|12p and (fT3")) and assuming 
that the site orbitals ipi = ^q(r — Hi) correspond to 
hydrogenic wavefunctions with quantum numbers q, the 
parameters of the Hamiltonian (|11[) can be expressed in 
terms of two-center molecular integrals as 

^ = ^q,'q _ 4<Sq, q X q ' q ) 

t = (£? q +io/-R)S , q , q -2Jf q , q 



x = - s q , q TO 



x qq 



_({/q^ . 



"TO], 



(15) 



where K, 



q,q 



A"^ q = Kq „, U is the on-site interaction, 



V and W are off-diagonal repulsion terms, and X is some- 
times called the density-dependent hopping or enhanced 
hopping rate [39[. The absolute value of the ratios be- 
tween different interaction parameters for atoms in np Zl 
ns and nd z 2 states {n = 30) are shown in Fig. O where 
we have plotted the parameters in Eq. (|15p obtained with 
Si = 0, that is, with Wannier functions represented by 
bare atomic orbitals. The orthogonalization procedure 
(i.e. replacing ipi by (pi) mainly affects the smallest pa- 
rameters X and W in the regions where the value of Sf 
is not small enough for the Wannier functions to be ef- 
fectively orthogonal. 

The interaction parameters calculated above corre- 
spond to the intermediate screening case, that is U > 

V > W and X ~ kJ, where k is a constant [13, El- 
The ratios between the parameters are similar to those 
found in realistic systems, e.g. between p-electrons in con- 
jugated polymers [40]. For different values of the angular 
momentum, these ratios differ mainly in the way they be- 
have at large distances. The hopping rates t are shown 
in Fig.[6l and, for the quantum numbers considered, vary 
between 10 2 and 10~ 4 GHz for intersite distances between 
150 and 250 nm respectively. 



(13) 

and the only two-electron interaction terms taken into 
account are U = V(i,i,i,i), X = V(i + 1, i, i, £), V = 
V(i, i + + and W = V(i, i + l,i,i + l) 

Electronic orbitals centered at different lattice sites are 
often considered to be orthogonal and equated to the 
Wannier functions of the crystal [37]. However, a better 
definition of the Wannier function centered at site i in 



IV. CONCLUSION 

We showed that optical lattices forming arrays of 
double-well potentials may be exploited to selectively 
photoassociate pairs of atoms to molecular states with 
binding energies of the order of 10 3 GHz, far larger than 
those of long-range molecules stabilized by dipole-dipole 
forces. These molecular states are expected to have equi- 
librium distances of the order of the typical lattice spac- 
ings and lifetimes several orders of magnitude larger than 



8 




R/a 



FIG. 6: (color online) Hopping rate in the z-direction between 
sites separated by a distance 7? for electrons in ns, np z , nd z 2 
states (n = 30). 



the timescales required to interrogate and characterize 
them by way of short laser pulses. 

We considered the possibility of collectively exciting 
ground-state atoms trapped in a regular lattice to a given 
Rydberg level such that the charge-density distributions 
of atoms located in neighboring sites overlap. Assuming 
that such a system could be realized, we calculated the 
typical interaction parameters between the Rydberg elec- 
trons and the hopping rate between sites. With a tem- 
perature well below the Fermi temperature and tunable 
interaction parameters, such systems might offer an in- 
teresting alternative approach to the simulation of Fermi 
systems. For instance, the ability to change the lattice 
spacing such that v=U/V>^otv<^ would offer the 
possibility of engineering a charge-density wave (CDW) 
or a spin-density wave (SDW), respectively. The prepa- 
ration of a SDW state could even be facilitated by using 
spin-changing collisions between ultracold 87 Rb ground- 
state atoms to produce a Neel-like state | Tlti ■ ■ ■) that 
has a spin arrangement identical to that of the SDW 
ground state along one spatial direction |42j. If realized, 
such a system would allow the observation of a phase 
transition between the SDW and CDW phases, by, for 
example, the measurement of the zero-frequency SDW 
susceptibility [53[ or the CDW structure factor both of 
which diverge linearly in their respective phase [43J . Fur- 
ther, the versatility of optical lattice setups may allow 
the excitation of Rydberg atoms with a given angular 
momentum, which potentially enables exotic quantum 
phase transitions to be engineered: indeed, the value of 
the electronic angular momentum not only influences the 
relative values of the system interaction parameters, but 
also in some cases their signs (see e.g. Refs. 0, EE 0] ) ■ 
This feature is particularly interesting, as it is the mech- 
anism behind e.g. the emergence of non-trivial and rich 
phase diagrams in doped cuprates (see e.g. Ref. [13, El])- 
Finally, if in addition the electron density could be con- 
trolled when producing one of these states, it would pro- 
vide the exciting possibility of turning the state of the 



electrons into a superconductor, as happens in the case of 
doping- induced superconductivity in cuprates [49]. The 
parameters calculated in this work apply to models where 
the dynamics of the electrons is restricted to one mode 
per site. 

In future works, it would be interesting to deter- 
mine the conditions for this assumption to be valid, and 
whether these conditions can be engineered with current 
technology. Results in this direction would allow excita- 
tion schemes to be devised, and also the determination 
of the control available to tune the density of electrons 
in the lattice. Also, dynamical calculations aiming at 
characterizing the lifetime of the lattice configuration af- 
ter the atoms have been excited to Rydberg levels would 
allow the evaluation of the time-scale available to inter- 
rogate the system. The exact calculation of this lifetime 
is an open problem whose solution is beyond the scope of 
the present paper. However, the dynamics of our system 
happens on a time scale of a few hundreds of picoseconds. 
Techniques for probing electron dynamics on the femto- 
second time scale exist in condensed matter systems (see 
e.g. Ref. [13|)). These methods might form the basis for 
resolving the motion of electrons in a Rydberg gas on 
pico- to femto-second timescales in future experiments. 
The use of our scheme for the purpose of quantum sim- 
ulation requires the lifetime of the lattice configuration 
to be longer than the typical timescales of the dynam- 
ics. As mentioned by Mourachko et al. in Ref. [5l|, the 
interatomic spacing (~ /im) between Rydberg atoms in 
a frozen gas varies very little (~ 3%) over a period of 
time of the order of 1 /is, some three orders of magnitude 
larger than the typical hopping times present in the sys- 
tem we have prop osed. Also, as already mentioned by Li 
et al. in Ref. [521 ] fixing the initial positions of the atoms 
will help reducing the motion of the ion cores. We there- 
fore believe that Rydberg gases created from ultracold 
atoms in an optical lattice provide a promising route to- 
wards the direct quantum simulation of interacting fermi 
systems. 
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APPENDIX A: REPLACEMENT RULES FOR 
THE EVALUATION OF TWO-CENTER 
MOLECULAR INTEGRALS WITH m = 0. 

Here we outline the methods we have used to solve 
molecular integrals using symbolic replacements. We 
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also provide all the necessary relations to implement the 
method for m = 0. 

The first step towards the evaluation of these integrals 
is to express the Coulomb potential in terms of a series 
of Legendre functions as 



1 



OO K 

^ ^ (/,- 

k=0 m=— k 



(k-\m\)l r(a) k 
fc+i 



i\)\r(b) 



P 1 ™ 1 {cos e^P 1 ™ 1 (cos 9 2 )e 



im(<p 1 -<t> 2 ) 



(Al) 



where r 12 is the distance between two points with spher- 
ical coordinates (ri,9i,cf>i), P^\x) are the associated 
Legendre functions [we use the notation P° (x) = Pk {x)\ , 
and r(a), r(b) are the smaller and larger of the quantities 
r\ and r 2 ; or using the von Neumann expansion 



_ oo k 



'"12 



k— m— — k 



(k - |m|)! 



1 2 



(fc + |m|)!_ 



(A2) 



where in this case the two points are expressed in ellip- 
tical coordinates (Aj, fii,(f>i), are associated Leg- 
endre functions of the second kind, and A (a) is the lesser 
and X(b) the greater of Ai and A 2 (see e.g. Ref. [Hj]). 

Since the integral e Il "^d0 vanishes if v is an integer 
different from zero, setting the quantum number m = 
considerably simplifies the evaluation of the integrals (j4|) 
and ([5]) using the relations (|A1|) and (|A2[) . 

Every two-center molecular integral apart from 
can be evaluated using Eq. (|A1[) . The angular part of 
these integrals will be non-zero for only a few terms; for 
instance, for q = q' = (n, 0,0), the angular part is non- 
zero for k = 0, and for q = q' = (n, 1, 0) for k — 0, 2. ft 
may be simplified by using the formula for the product 
of surface harmonics: 



Pt x (cos 9)P? 2 (cos ( 



[Cj lA ] 2 P,(cos0), (A3) 

j=\h-e 2 \ 



where C- 1 



= C(lx,t2, j;0,0,0) are Clebsch-Gordan co- 
efficients [31| . The integration of the radial part over the 
coordinates of the first electron necessitates the evalua- 



tion of integrals of the form 



v(ci) k 



(A4) 



where a q (?") is the radial part of the wavefunction ([3]), and 
ri is the coordinate of the i th electron. The functionals 
J q t q / jfc (r 2 ) are defined as 



^ kfa) = (lA2 +i ) / dn rir^r^in) (A5) 
Jo 

and 



lL>,k(r 2 ) = r« dnr^r^MnKKn). (A6) 



For fc 
rules 



both (|A5|) and (|A6|) can be solved using the 



R3 : / dre~ ar r 



ar^k 



R4 : / dr e~ ar r k 



[k\-T k+l (ar)}/a k+1 
T k+ i(ar)/a k+1 (A7) 



If the sum of Eq. (|A1[) contains only fc = terms, the 
symbolic form of the integral is obtained by replacing 
the spherical coordinates of the remaining electron by el- 
liptical coordinates and using the replacement rules of 
Eqs. (|A7[) and Table [IJ In integrals requiring terms as- 
sociated with k > in the sum (IA1[) . the functionals 
associated with k = dominate [361] , For higher values 
of k, the functional (|A6|) can always be solved, and as 
an approximation we have dropped the terms in (|A5|) 
that could not be solved analytically using the integra- 
tion rules mentioned above. 

Because it involves associated Legendre functions of 
the second kind, the evaluation of W^ ^, is more prob- 
lematic. We used the Rodrigues formula to expand these 
functions in terms of sums of Legendre polynomials and 
logarithms, and then solved the polynomial part by ap- 
plying the replacement rules mentioned above. Solv- 
ing the logarithmic part requires the integration over 
log(l ± x) for the first set of coordinates, and then expo- 
nential integral functions for the second set. 



[1] J. M. Hutson and P. Soldan, Int. Rev. in Phys. Chem. 

25, 497 (2006). 
[2] G. Thalhammer, K. Winkler, F. Lang, S. Schmid, R. 

Grimm, J. H. Denschlag, Phys. Rev. Lett. 96, 050402 

(2006) . 

[3] B. Knuffman and G. Raithel, Phys. Rev. A 75, 053401 

(2007) . 

[4] D. DeMille, Phys. Rev. Lett. 88, 067901 (2002). 



[5] C. H. Greene, A. S. Dickinson, and H. R. Sadeghpour, 

Phys. Rev. Lett. 85, 2458 (2000). 
[6] A. A. Khuskivadze, M. I. Chibisov, and I. I. Fabrikant, 

Phys. Rev. A 66, 042709 (2002). 
[7] S. M. Farooqi et al, Phys. Rev. Lett. 91, 183002 (2003). 
[8] C. Boisseau, I. Simbotin, and R. Cote, Phys. Rev. Lett. 

88, 133004 (2002). 

[9] J. Stanojevic et al, (2008L larXiv:080L2386V l. 



10 



T. Gericke et al, (2008), (Preprint: larXiv:0804.4788) . [33 
R. Saers et al, The European Physical Journal Applied 
Physics 42, 269 (2008). [34 
K. I. Petsas, A. B. Coates, and G. Grynberg, Phys. Rev. 
A 50, 5173 (1994). [35 
M. Lewenstein et al., (Preprint cond-mat/0606771 ) 
(2006). ' [36 

M. Kohl and T. Esslinger, Europhysics News 37, 18 
(2006). [37 
Z. Idziaszek, L. Santos, and M. Lewenstein, Phys. Rev. 
A 64, 051402(R) (2001). [38 
T. Cubel et al, Phys. Rev. A 72, 023405 (2005). 
J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V. [39 
Porto, Phys. Rev. A 73, 033605 (2006). 

P. W. Atkins and R. S. Friedman, Molecular quantum [40 
mechanics (Oxford University Press, Oxford, 2005). 
J. C. Slater, Quantum Theory of Molecules and Solids [41 
(Mc-Graw-Hill, New- York, 1974), Vol. 1. 

J. Stanojevic et al, Eur. Phys. J. D 40, 3 (2006). [42 

C. A. Coulson, Valence (Oxford University Press, Oxford, [43 
1952). [44 
W. Li, I. Mourachko, M. W. Noel, and T. F. Gallagher, 
Phys. Rev. A 67, 052502 (2003). [45 

A. V. Stolyarov and M. S. Child, J. Phy. B: Atomic, 
Molecular and Optical Physics 32, 527 (1999). [46 

B. G. Wilson, C. Iglesias, and F. Rogers, Phys. Rev. A [47 
38, 4633 (1988). [48 
M. Marinescu and A. Dalgarno, Phys. Rev. A 52, 311 
(1995). [49 
M. Marinescu, Phys. Rev. A 56, 4764 (1997). [50 
R. J. Buehler and J. O. Hirschfelder, Phys. Rev. 83, 628 [51 
(1951). [52 
M. P. Barnett, International J. Quantum Chem. 95, 791 [53 
(2003). 

J. Fernandez Rico, R. Lopez, G. Ramirez, and C. 

Tablero, Phys. Rev. A 49, 3381 (1994). 

M. P. Barnett, J. Symb. Comput. 42, 265 (2007). 

M. P. Barnett, J. Chem. Phys. 113, 9419 (2000). 

D. Jaksch et al, Phys. Rev. Lett. 81, 3108 (1998). 



T. F. Gallagher, Rydberg atoms (Cambridge University 
Press, Cambridge, 1994). 

E. Toth, A. M. Rey, and P. B. Blakie, (2008), (Preprint: 
arXiv:08032922). 

M. Kohn and T. Esslinger, Europhysics News 37, 18 
(2005). 

J. Appel, M. Grodzicki, and F. Paulsen, Phys. Rev. B 
47, 2812 (1993). 

N. W. Ashcroft and N. D. Mermin, Solid State Physics, 
(Brooks/Cole Thomson Learning, Philadelphia, 1976). 
A. Painelli and A. Girlando, Phys. Rev. B 39, 2830 

(1989) . 

F. Marsiglio and J. E. Hirsch, Phys. Rev. B 41, 6435 

(1990) . 

J. T. Gammel and D. K. Campbell, Phys. Rev. Lett. 60, 
71 (1988). 

D. Baeriswyl, P. Horsch, and K. Maki, Phys. Rev. Lett. 
60, 70 (1988). 

S. Trotzky et al., Science 319, 295 (2008). 

A. M. Rey et al., Phys. Rev. A 72, 023407 (2005). 

E. Kaxiras, Atomic and electronic structure of solids 
(Cambridge University Press, Cambridge, UK, 2002). 

I. S. Elfimov, S. Yunoki, and G. A. Sawatzky, Phys. Rev. 

Lett. 89, 216403 (2002). 

J. E. Hirsch, Phys. Rev. B 48, 9815 (1993). 

D. Duffy and A. Moreo, Phys. Rev. B 52, 15607 (1995). 

M. Raczkowski, R. Fresard, and A. M. Oles, Europhysics 

Letters 76, 128 (2006). 

S. Sachdev and S.-C. Zhang, Science 295, 452 (2002). 
A. Cavalleri et al., Phys. Rev. Lett. 95, 067405 (2005). 
Mourachko et al, Phys. Rev. A 70, 031401(R) (2004) 
Li et al., Phys. Rev. Lett. 94, 173001 (2004) 
The SDW susceptibility and the structure 
factor are defined respectively by x(q) = 
■ 1/T dTE w Ki(r)/Sj(0) and S(q) 
l i( R i- R j){ ni nj), where Ki(t) = n^(i) 
and m. a = cj„ Ci 



nil 



